Self organized criticality in a sandpile model with threshold dissipation 
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We study a nonconservative sandpile model in one dimension, in which, if the height at any site 
exceeds a threshold value, the site topples by transferring one particle along each bond connecting 
it to its neighbours. Its height is then set to one, irrespective of the initial value. The model shows 
nontrivial critical behavior. We solve this model analytically in one dimension for all driving rates. 
We calculate all the two point correlation functions in this model, and find that the average local 

\f} • height decreases as inverse of the distance from the nearest boundary and the power spectrum of 

^\ ' fluctuations of the total mass varies as 1//. 
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Self-organized Criticality (SOC) was proposed by Bak, Tang and Wiesenfeld [1] to explain the widespread occurrence 
K-5 1 of fractal structures and 1/ f noise in nature. They suggested that extended dissipative systems evolve to a critical 
^| ! state by a self-organizing process and proposed the sandpile model as a prototype of SOC systems. Since then several 
sandpile-like stochastic cellular automata models with threshold dynamics have been studied to understand the origin 
of criticality. In models with bulk mass conservation [2-5], the balance of rate of driving and the rate of loss of 
t-H , particles through the boundary leads to a divergence of the average avalanche size in the steady state. However, the 
origin of criticality in nonconservative sandpile models is still not well understood [6,7]. In particular, the nature 
of nonconservation which can preserve the criticality of the model are not known. Analytic proofs of criticality in 
nonconservative sandpile models are lacking even in one dimension. 

In this paper we study a discrete nonconservative sandpile model [9] in which mass dissipation (loss of a particles) 
occurs if the height at a site exceeds a threshold value. In the steady state, our model organizes such that mass 
dissipation is minimized. As a result, it shows nontrivial critical behavior. We calculate its spatial and temporal 
ON • correlation functions analytically in one dimension. The average height profile in the steady state of this model has a 
nontrivial power law dependence on the distance from the nearest boundary. We find that the power spectrum of the 
fluctuations of the total mass of the sandpile shows 1// behaviour. These power laws are robust and do not depend 
on the details of the model. The model is critical even for finite driving rate (unlike the forest fire model [8] , which 
is critical only for infinitesimal driving rate). To our knowledge, this is the first nonconservative model of SOC which 
has been shown to be critical for all driving rates. 

The model is defined as follows: There is an integer variable hi (the number of particles), at each site i. The 
particles are added to the system at randomly selected sites. If hi exceeds the threshold height h\, then hi is set to 1 
and one particle is transferred to each neighbor of i. This process occurs at all sites in parallel. We choose h\ to be 
• i-H . equal to the coordination number for i inside the bulk and to be greater than the coordination number for i on the 
boundary. Note that this model is not abelian [2]. In the Abelian sandpile model (ASM), hi decreases by hi after the 
toppling, whereas in this model hi is set to 1 irrespective of its initial value. 

The dynamics of this model resembles the dynamics of the continuous stick-slip model discussed by Feder and Feder 
(FF) [6]. The FF model has a real variable Ui at each site i, which grows slowly with time. If Ui exceeds U c , then Ui 
is set to zero and Uj is increased by 1 for all neighbors j of site i. Let us define hi as the integer part of (ui + 1) and 
hi as (U c + 1). Then clearly, the updating rule can be written in terms of hi and it is identical to the toppling rule 
of our model. However, the external driving force in FF model is deterministic and not noisy as in our model. Note 
that the argument given above crucially depends on the fact that the transfer to the neighbours does not depend on 
the local variable. Therefore, it does not hold for the earth-quake model [7], where the transfer is proportional to the 
local variable. 

To illustrate the difference between this model and the ASM, consider the model on a square lattice and choose an 
initial configuration in which the four corners of a given plaquette have all heights equal to 4 (threshold height) . Now 
add a particle at the lower left corner of the plaquette, say at («, j). First (j, j) topples then (i + 1, j) and (i, j + 1) 
topple. After this the height at (i + + 1) becomes 6. So far there is no dissipation (just as in the ASM). When 
(i + 1, j + 1) topples, unlike in the ASM, 1 particle is lost. Note that this model is the same as the ASM if the 
underlying graph does not have a loop. 

In spite of the non-abelian character, the set recurrent configurations of this model Si is the same as that of the 
Abelian sandpile model S2. To prove this, we first show S2 C Si. Let us add a particle in Co € S2 and let it relax 
by the rules of the ASM. Let the final configuration be C\. Now we add a particle at the same site in Co and let 
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it relax by the rules of this model. If there is no dissipation then the final configuration is the same as C\ . If there 
is dissipation, then we compensate for the dissipation by adding particles one by one at sites where dissipation has 
occurred and letting the system relax. This process is repeated until no site is left where dissipation has occurred. 
Since in the ASM the final configuration does not depend on the order of the topplings, this sequence of additions 
and relaxations in this model must lead to C\. Thus, in this model, there is a finite probability of transitions from 
one recurrent configuration of the ASM to the another. This implies that £2 C S±. One can easily check using 
the argument of Ref. [2] that the forbidden configurations of the ASM are also forbidden in this model. Therefore, 
Si = £2. However, unlike the abelian case, the probabilities of occurrence of different recurrent configurations in the 
steady state of this model need not be equal. 

Let us analyse this model in one dimension. The model on a simple linear chain is not different from the ASM, 
because this lattice has no loops. Hence, we consider this model on decorated one dimensional chains, formed by 
joining unit cells which have loops (see Fig. 1). We have earlier studied ASM's on this type of chain and found 
interesting finite-size scaling behavior [10], which differs from that seen in simple linear chains [11]. 

Consider first the chain of doublets (case A, Fig. 1). We label the doublets by integers i — — I to I and the sites 
inside the doublet by j = 1 to 2. The size of the chain L = 21 + 1. We generalize the rule such that after the toppling 
one particle is transferred along each bond connecting the site to its neighbors, and take the threshold height equal 
to 3 for all sites. The recurrent configurations of this model are characterized using the burning algorithm as in the 
ASM This algorithm is specified by the following rule: A site is burnt if its height is greater than the number of bonds 
joining it to its unburnt neighbours (see [2], for detail). A stable configuration is recurrent if and only if all the sites 
are eventually burnt. 

In this algorithm the sites can be burnt in any order. We let the burning start from the left boundary and hold right 
boundary unburnt. The point at which burning from this boundary stops is called the break point (BP). Afterwards, 
the right boundary is burnt and subsequently the remaining sites. Thus the allowed values of (hnjhiz) for i to the 
left of the BP are (3, 3) and (3, 2), and those for i to the right of the BP are (3, 3) and (2, 3). The allowed values at 
the BP are (1,3), (3,1) and (2,3). 

In any side the avalanche spreads to the first doublet which cannot be burnt from that side. For example, if the 
particle is added to the left of the BP the avalanche spreads to the BP on the right and to the first doublet of type 
(3, 2) on the left. For a typical avalanche the distance between the site at which the particle is added and the BP is 
O(L). Therefore, the spread of avalanche is O(L). 

When the avalanche crosses a doublet of type (3,3), one particle is dissipated. Thus, each avalanche wipes out 
(3, 3) type of doublets from a region of order L, leaving a small density of the (3, 3) type doublets. This minimizes 
the possibility of dissipation. Also, since the steady state configurations are dominated by (3, 2) type doublets ((2, 3) 
type doublets) on the left (right) of the BP, the leftward (rightward) spread of the avalanche starting from the left 
(right) of the BP is small (see Fig. 2). 

The BP shows an interesting stochastic dynamics. After the avalanche, it moves towards the starting point of the 
avalanche, by a distance of order 1. If i is the position of the BP then the probability that an avalanche starts from 
the left of the BP is (I + i)/(2L) to leading order (ignoring the density of (3, 3) type of doublets), and that it starts 
from the right of the BP is [l — i)/(2L). Hence, the mean displacement of the BP after the avalanche scales as i/(2L). 
The fluctuation about this mean value is of order 1. The dynamics of the BP can thus be described by an equation 
of the following form 

dx x 

where x is the scaled variable i/L, and dt is of the order of the time interval at which avalanches hit the break point. 
The noise r)(t) <~ 1/L, and is (5-correlated. From the above equation it follows that the number of avalanches required 
to reach the steady state is of order L, and the asymptotic distribution of the BP goes as exp(— cx 2 L), where c is a 
constant. The width of this distribution can be ignored in the large L limit. A typical avalanche then starts from the 
point at which particle is added (source point) and terminates at the center. If the source point is i then the linear 
extension s and the duration t of an avalanche are equal to \i — L/2\. Averaging over i we get 

Prob(A) - 2/L for X < L/2 , (2) 

and otherwise, where X = t, s. 

Now we calculate the correlation functions in this model. The model has two intrinsic time scales, the time taken 
for one toppling (defined as one time step), and the time interval T between two consecutive addition of particles. 
Let us first consider the case, T > L, such that there is no overlap between two consecutive avalanches (see Eq. (2)). 

It is convenient to use a two state variable Sj , because there are only two allowed configurations of each doublet 
(except the break point). For the (3,3) type doublet, Sj = 1, otherwise Sj = 0. Consider the evolution Sj for 
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VL "C i <C l- We can assume that the BP is on left of i. If the avalanche starts from the right of i, then it crosses i 
to reach the BP, setting Si — 0. The probability of this transition goes as (I — i)/(2L). The transition of Si from to 
1 occurs, if the particle is added at the left site of the i-th doublet, or if it is added to the left of the i-th doublet, but 
the avalanche reaches the (i — l)-th doublet so that the left site of i-th doublet receives a particle. The probability 
of this transition goes as 1/L. The corresponding probabilities for -I < i « —VL, can be obtained by replacing i 
by (— i) and I by (— I). Thus the probability P n (i), that Si = 1 after additions of n particles, satisfies the following 
equation (to leading order) 

Pn+l(i) = \ [1 - Pn(i}\ ~ (l - ^) Pn{i) , (3) 

where r{i) is the distance of i from the nearest boundary. A straightforward calculation using the above equation 
gives 

(Si)~2/r(i) , (4) 

and the autocorrelation 

(sdto) S dt +t))~-±- ) (l-T- 1 ) t/T , (5) 

for t > T, where n = and (. . .) denotes average over to- The value of Si at time t is denoted by Si(t). Note that 
the correction to Eq.(3) coming from the nonzero density of Si (Eq. (4)) vanishes in the large L limit. The steady 
state of this model is not translationally invariant, but is self similar in space, i.e. the density and the relaxation time 
Ti depend on its distance from the nearest boundary as a power law. In Fig. 3, we plot the density of Si against the 
distance from the left boundary. Note the excellent agreement of our numerical data with Eq. (4) . 

To calculate the correlation between Si(t) and Sj(t'), consider the evolution of the joint probability of Si and Sj. 
There are four possible values of (si,Sj). Let the four-vector, P n (i,j) 7 denote the joint probability of (si,Sj) after 
addition of n particles. Its evolution can be described by the following linear equation (to leading order) 
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(i,j) = G(i,j)P n (i,j) , (6) 



where G(i,j) is the (4 x 4) matrix giving the transition probabilities of Si and Sj. This can be calculated in exactly 
the same way as the transition probability of Si calculated above. For example, (1,1) goes to (0,0) if the avalanche 
passes over both the doublets i and j. If \f~L < i < j < i then the probability of this transition goes as r(j)/(2L). 
Other transition probabilities can also be calculated in the same way. From Eq. (6) one can find the steady state 
correlation functions. For i and j on the same side of the center of the chain 

(«i(*0) Sj(t + t)) c ~ ^ (1 - T7 l ) t/T , (g) 

for r(i) > r(j) , 

(«i(to)gj(*o+*))c~ - Tmruj^m ( 1 ~ r r 1 ) t/T . (9) 

for r(i) < r(j). 

The subscript c refers to the connected part of the correlation function. If the i and j are on different sides of the 
center then the correlation can be ignored. Note that, in Eq. (8), the first factor gives the equal time correlation, 
and the second factor gives the relaxation of Sj. In Eq. (9) the s,(£o) and Sj(to + 1) are anti-correlated because an 
avalanche which makes Sj = 1 sweeps over site j, setting Sj =0. 

The autocorrelation of mass (sum of heights) can be obtained using Eq. (5), (8) and (9). The power spectrum of 
mass fluctuation S(f), is obtained by taking the real part of the Fourier transform of the autocorrelation 

1 2tt 

S(f)~j , for — «/«2tt . (10) 
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In Fig. 4, we show the log-log plot of the numerical value of S(f) versus /. This is obtained by taking the square of 
the Fourier transform of the time sequence of length 5000, for a chain of size 5000. We find that at low frequencies 
the data shows 1/f behaviour. The scaling region increases with the increase in the lattice size. Since we have taken 
only the slowest relaxation mode in Eq. (8), (9) and (5), the agreement is not good in the high frequency regime. 

Consider now the case in which particles are added at a faster rate, i.e., 1 < T « L. We note that the interaction 
between the avalanches can be ignored because the avalanches propagate with the same velocity towards the BP (see 
Fig. 2). Thus the steady state properties of the model as described by Eq. (7) and (4) and the autocorrelation given 
in Eq. (5) remain unchanged in this case. Since the avalanche takes 2\j — i\ time steps to travel from i to j, there is 
no anti-correlation between Sj(£o) and Sj(to + 1) for t < 2\j — i\ (see Eq. (9)). However, the contribution to the mass 
autocorrelation coming from the this term can be ignored for \i — j\ > \/L, or for t > \[L. Thus for / < 2ttT/^/L, 
the power spectrum S(f) ~ 1/f. 

To check the robustness of our results we have studied the model on the diamond chain (case B in Fig. 1). An 
avalanche in this case has same structure as in case A, i.e., it spreads to the BP on one side, to a distance of 0(1) 
on other side and the BP is confined in a region of 0(vL). As a result, the power laws of the average height and the 
power spectrum of mass fluctuation are also the same. 

The two special doublet configurations in this model, which stop the avalanches, are reminiscent of the 'trough' 
and 'trap' of the 'singular diffusion' type model studied by Carlson et al [5]. In that case too there is a power law 
decay of the height profile which follows from the singularity of the effective diffusion coefficient of the particles as a 
function of the coarse grained density [4,12]. However, in our model the mechanism of self-organization is different as 
the diffusion coefficient does not diverge. 

To summarize, we have determined exactly the critical behavior of a non-abelian nonconservative sandpile model 
in the one dimension. The critical steady state shows spatial structures. We determine time dependent correlation 
functions for finite driving rate, and show that fluctuations of the mass of the sandpile have a 1// spectrum. 

I thank Prof. Deepak Dhar for his useful suggestions and in particular, for pointing out that the set of the recurrent 
configurations of this models is same as that of the ASM. I thank Gautam I. Menon and P. Lakdawala for a critical 
reading of the manuscript. 
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Captions 

Fig. 1: The one dimensional chains formed by joining (A) doublets, (B) diamonds and (C) single sites. 
Fig. 2: The evolution the sandpile. The dot denotes an toppling event. 

Fig. 3: The log-log plot of the mean height of a doublet (minimum height is subtracted) versus the distance from 
the nearest boundary r. The solid line shows the analytic expression. 

Fig. 4: The log-log plot of power spectrum of total mass fluctuation S(f) vs frequency /. For comparison the 
theoretical curve is shown as the solid line. 
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